Donor and acceptor levels of organic photovoltaic compounds from first principles 
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Accurate and efficient approaches to predict the optical properties of organic semiconducting 
compounds could accelerate the search for efficient organic photovoltaic materials. Nevertheless, 
predicting the optical properties of organic semiconductors has been plagued by the inaccuracy or 
computational cost of conventional first-principles calculations. In this work, we demonstrate that 
orbital-dependent density-functional theory based upon Koopmans' condition [Phys. Rev. B 82, 
115121 (2010)] is apt at describing donor and acceptor levels for a wide variety of organic molecules, 
clusters, and oligomers within a few tenths of an electron-volt relative to experiment, which is 
comparable to the predictive performance of many-body perturbation theory methods at a fraction 
of the computational cost. 



I. INTRODUCTION 

Despite the positive attributes of organic materials for 
next-generation photovoltaics, the effective deployment 
of organic photovoltaic (OPV) cells poses fundamental 
problems at the molecular scale. Besides enhancing the 
durability of organic compounds, the central requisite to 
the mass-market viability of OPV modules is to raise 
their power conversion efficiency (PCE) beyond 10-15%^^ 
Comparatively, the efficiency of current OPV architec- 
tures (Fig.[T]) barely exceeds 9-10% PCE under standard 
test conditions.^ Nevertheless, the margin for improve- 
ment is vast due to the exceptional chemical versatility 
of organic materials,^ and accurate approaches to predict 
the charge-transfer and optical properties of new com- 
pounds could allow for extensive screening of promising 
semiconducting organics.^ 

However, conventional density- functional theory ap- 
proximations in their static (DFT) and time-dependent 
(TDDFT) formulations^ suffer from severe limitations 
in predicting the electronic structure of semiconduct- 
ing materials,^ precluding the quantitative and even 
qualitative elucidation of relevant photovoltaic mecha- 
nisms. Notably, conventional DFT and TDDFT ap- 
proximations tend to systematically destabilize occupied 
states and overstabilize unoccupied states in organic and 
organo metallic complexes, leading in particular to the in- 
correct description of the electronic properties (charged 
excitations) and optical properties (neutral excitations) 
of donor-acceptor dyads, extended polymer molecules, 
and heavy-metal sensitizing dyes.^^^ To overcome DFT 
Hmitations, one privileged route has been to resort to 
many-body approaches, namely, to many-body pertur- 
bation theory approximations such as GW .^^i^^ The GW 
method, while accurate in predicting electronic spectra. 
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FIG. 1: Schematics of an organic solar cell showing a bilayered 
donor-acceptor heteroj unction (top). Bilayered, nanopat- 
terned, and mixed-phase heteroj unctions with typical power 
conversion efficiencies (bottom). 



is much more expensive than DFT although considerable 
progress has been achieved in reducing the cost of GW 
calculations 

In parallel to many-body perturbation theory, less 
expensive orbital-dependent density-functional theo- 
ries (OD-DFTs) represent promising alternatives At 
present, the most widely used OD-DFT methods are hy- 
brid density- functional theory approximations.^^ Hybrid 
functionals in their linear- admixture or range-separated 
forms have been shown to improve upon conventional 
DFT in predicting electronic properties. ^^^^^^^ Alterna- 
tively, self-interaction corrections to density-functional 
theory approximations^^ represent a second category of 
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OD-DFT methods that are now rapidly growing in recog- 
nition due to their improved accuracy and lower compu- 
tational cost. In self-interaction-corrected approaches, 
the total energy of the system is expressed explicitly in 
terms of individual orbital densities to rectify nonphysi- 
cal errors inherent in orbital-independent DFT. Success- 
ful applications of advanced self-interaction corrections 
have appeared with an accuracy potentially approach- 
ing that of many-body methods in predicting electronic 
levels and electrical responses^ — 

In this work, we highlight the predictive accuracy of 
the recently developed self-interaction correction based 
upon Koopmans' theorem (Koopmans-compliant OD- 
DFT)^i2^ in capturing the electronic levels of OPV 
compounds. Electronic levels represent reliable PCE 
indicators^ and have successfully served as fundamen- 
tal inputs for first-principles combinatorial screening 
of organic semiconductors.^'^ We demonstrate that 
Koopmans-compliant OD-DFT is apt at describing donor 
levels within 0.1-0.4 eV and acceptor levels within 0.2- 
0.6 eV relative to experiment, which is comparable to 
the precision of the GW method. In its simplest formu- 
lation, the method can be trivially implemented in con- 
ventional electronic-structure codes with an algorithmic 
cost potentially lower than that of hybrid-DFT function- 
als as it does not require evaluating exchange terms. Fur- 
thermore, the method holds promise for accurate time- 
dependent extensions based upon the improved descrip- 
tion of the underlying electronic spectrum. 

This work is organized as follows. First, we underscore 
the significance of Koopmans' theorem in understand- 
ing the charge-transfer behavior of donor-acceptor com- 
plexes. Second, we outline existing electronic-structure 
approaches that aim at enforcing Koopmans' theorem. 
Third, we present the Koopmans-compliant OD-DFT 
method to correct the nonphysical tendency of conven- 
tional DFT approximations to destabilize occupied states 
and overstabilize unoccupied states. Last, we demon- 
strate the efficiency and accuracy of the method for fam- 
ilies of organic compounds of interest to photovoltaics. 



II. METHOD 

A. Koopmans' theorem 

Koopmans-compliant functionals aim at imposing 
Koopmans' theorem in DFT approximations. Koop- 
mans' theorem enables one to equate orbital energies 
(that is, the expectation values of the effective Hamil- 
tonian) with total energy differences, which correspond 
to withdrawing an electron from a stationary electron 
state. 

The relevance of Koopmans' theorem on the accu- 
racy of effective-potential theories in predicting elec- 
tronic spectra has been recognized in several theoretical 
and computational studies j^ i^^i^^ — In addition, Koop- 
mans' theorem is the central condition to correctly de- 
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FIG. 2: Energy of a positively charged pentacene-fullerene 
dyad in the infinite separation limit (black) and energies of 
the isolated fractionally charged pentacene (light blue) 
and fullerene Cg^^ molecules (dark blue) as a function of 
the transferred charge q within the LSDA and Koopmans- 
compliant descriptions. (Energies are given relative to the 
fully ionized individual molecules.) 



scribe electron transfer within donor-acceptor pairs. To 
illustrate this fact, we consider a positively charged 
pentacene-fullerene dyad in the limit of infinite inter- 
molecular separation where electronic interactions be- 
tween the two molecules can be neglected. 

Figure [2] depicts the energies of the isolated pentacene 
oligomer and fullerene cluster. The total energy of the 
infinitely separated molecular pair as a function of the 
transferred charge q is also reported. Within the local 
spin-density approximation (LSDA) [Fig.[2fa)], the three 
curves exhibit a marked nonlinear dependence on q^ re- 
flecting the fact that LSDA violates Koopmans' theorem. 
In other words, the derivative of the energy dE{yi^) /dq 
of an individual fractionally charged cation X^, which 
corresponds to the opposite energy of its highest occu- 
pied orbital, departs significantly from the finite energy 
difference EilC^) — £^(X), which corresponds to the ion- 
ization potential of the neutral molecule X. More pre- 
cisely, the dependence of the energies is found to be con- 
vex, causing the dyad to be most stable at a fractional 
charge ^min close to |. This observation is in quali- 
tative contradiction to the expectation that the charge 
should fully localize onto the strongest electron accep- 
tor (namely, fullerene). In contrast to this nonphys- 
ical behavior, Koopmans-compliant energy curves are 
linear [Fig. [2fb)] as a result of Koopmans' condition 
d^(X^)/d(7 = ^(X+) - EiX), leading to the expected 
stabilization of the fully transferred charge (g'min = !)• 

Similar observations would be made for other donor- 
acceptor complexes. Moreover, at finite separation where 
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electronic coupling between the two molecules must be 
taken into account, LSD A would also significantly under- 
estimate the transferred charge. These results highlight 
the central importance of imposing Koopmans' theorem 
on the individual subsystems to correctly describe charge 
transfer in donor-acceptor pairs. 



B. Existing methods 

At present, there exist different DFT-based methods 
that aim at imposing Koopmans' theorem j^^i^^i^^" — The 
self-interaction correction of Perdew and Zunger^^ can be 
regarded as one of the earliest approaches to restore the 
physical interpretation of orbital energies as ionization 
energies. In this approach, the Kohn-Sham density func- 
tional E^^ is augmented with corrections that remove 
the self-interaction of the electronic wave functions in the 
one-electron limit. In explicit terms, the Perdew- Zunger 
(PZ) correction to the Kohn-Sham DFT functional reads 
^pz ^ ^KS _ ^.^E^^^Ip,,], where pi, stands for the 
density of the one-electron state ipi^ and £^hxc denotes 
the Hartree plus exchange-correlation component of the 
Kohn-Sham functional. The PZ approach is more ex- 
pensive than local and semilocal DFT approximations 
but appreciably less costly than hybrid-DFT approxima- 
tions. It rectifies the tendency of local and semilocal DFT 
to destabilize occupied orbitals in atoms but overcorrects 
the error for molecules, leading to an overestimation of 
ionization energies of up to 1-2 eV. Also, the PZ correc- 
tion leaves the energies of empty levels unchanged, mean- 
ing that acceptor levels remain largely underestimated 
(that is, excessively negative). Furthermore, the lack of 
balance in the PZ correction of occupied and unoccupied 
orbitals leads to the underestimation of the total ener- 
gies of molecules and the underestimation of interatomic 
distances. 

A solution to the lack of balance in the PZ correc- 
tion has been proposed by Baer et al. in the form of a 
range-separated hybrid functional [the Baer-Neuhauser- 
Livshits (BNL) functional] in which long-range interac- 
tions are treated at the explicit-exchange level, whereas 
short-range interactions are described at the semilocal 
density-functional level^^i^ In order to impose Koop- 
mans' theorem within BNL, the original empirical pa- 
rameter 7 that determines the range of the electrostatic 
separation (1/r = erf(7r)/r + erfc(7r)/r) is optimally 
tuned to impose the agreement between orbital energies 
and ionization energies for the frontier orbitals of the 
system, ^^^^^ leading to considerable improvement in pre- 
dicting electronic and optical spectra. ^^^^^ In terms of 
computational cost, BNL calculations are as expensive 
as hybrid-DFT calculations. 

A recently proposed approach that targets mini- 
mal computational cost is the method of Lany and 
Zunger, ^^^^^ which consists of rectifying deviations from 
Koopmans' theorem through on-site projection terms 
that shift the electronic levels in order to impose the 



correspondence between orbital energies and ionization 
energies. Because it only relies on the calculation of on- 
site occupations, the Lany-Zunger method does not cause 
any noticeable increase in computational cost relative to 
the underlying DFT approximation. However, it requires 
to preselect atomic orbitals to construct the on-site pro- 
jectors, making its application less systematic than the 
PZ correction and optimally tuned BNL method. 



C. Present method 

The Koopmans-compliant method that we have intro- 
duced in Ref. 28 provides a promising alternative to ex- 
isting approaches. The method is systematic and allies 
conceptual simplicity, computational performance, and 
predictive accuracy, thereby permitting the precise and 
efficient electronic-structure description of a wide variety 
of compounds. 

The motivation for introducing the Koopmans- 
compliant method originates from the practical obser- 
vation that determining ionization potentials (IPs) of 
molecules as differences of ground-state energies In = 
En -I — En (where In and En denote the IP and to- 
tal energy of the system of TV electrons) yields predic- 
tions in very good agreement with experiment within 
conventional DFT approximations, namely, the LSDA 
approach and semilocal generalized-gradient approxima- 
tions (GGAs) — in the literature, the procedure that 
consists in evaluating ionization energies from differences 
of DFT total energies is commonly known as the ASCF 
method. For electron affinities (EAs), ASCF predic- 
tions An = En — En-\-i are also found to be in agreement 
with experiment with the central caveat that ASCF pre- 
dictions for EAs are only possible insofar as the system 
can accommodate the added electron within the approx- 
imation used; this condition on the stability of anionic 
states is necessary to calculate the ground-state energy 
of negatively charged systems. In practical terms, fail- 
ure to stabilize anions is frequent within approximate 
DFT^^ and arises from spurious self-interaction whereby 
an electron can interact with itself through the nonphys- 
ical contribution from its own charge density to the ef- 
fective Kohn-Sham Hamiltonian.^^ 

Therefore, in order to obtain reliable predictions for 
orbital levels (especially, for frontier donor and accep- 
tor levels that determine the charge-transfer properties 
of the system), it is crucial to satisfy the following con- 
ditions: (a) the correspondence between the effective or- 
bital levels and ASCF total-energy differences must be 
imposed; (b) the established accuracy of ASCF energies 
must be preserved; and (c) the calculations must circum- 
vent the inability of approximate functionals to stabilize 
negatively charged systems. Fulfilling these requirements 
is the object of the non-Koopmans OD-DFT correction 
that is described below. 

In presenting the non-Koopmans correction, it is im- 
portant to distinguish between electronically adiabatic 
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FIG. 3: Diagram of the non-Koopmans energy n2p, which 
quantifies departure from linearity for the adiabatic ionization 
curve of the 2p state of carbon within LSD A. 



(i.e., unrelaxed) ionization processes, whereby orbitals 
are kept frozen upon ionizing the system, and diabatic 
(i.e., relaxed) ionization processes, whereby orbitals are 
allowed to rearrange self-consistently. For adiabatic ion- 
ization, the correspondence between orbital energies and 
total energies is referred to as the restricted Koopmans' 
theorem, whereas for diabatic ionization, it is termed the 
generalized Koopmans' theorem. In what follows, we 
consider the adiabatic case before addressing the diabatic 
case. 

Bearing in mind this distinction, the starting point of 
the method consists in defining an adiabatic measure of 
the deviation from Koopmans' theorem in terms of an 
orbital-dependent energy functional. We thus consider 
for the moment an unrelaxed process that corresponds 
to depleting the stationary state i/ji^. For this process, 
the energy measure that monitors deviations from the 
restricted Koopmans' theorem can be written as the non- 
Koopmans energy 

nUf)= f^^' df'[eUf)-eUf% (1) 

^0 

first introduced by Perdew and Zunger in discussing 
the nonlinear behavior of the DFT total energy upon 
ionization. In Eq. ([!]), the superscript notation {-j^ in- 
dicates that the orbitals are kept unrelaxed along the 
ionization path, fia- denotes the initial occupation of the 
ionized orbital, and ef^ stands for the expectation value 
of the effective Kohn-Sham Hamiltonian for the ionized 
state ipia- along the unrelaxed curve. In other words, e^^ 
equals the derivative of the unrelaxed total energy Ef^ 
as a function of the occupation of the ionized state (the 
restricted Janak theorem). Alternatively, the unrelaxed 
non-Koopmans energy can be written in explicit terms 



as 

+ j dr^xca (r; [/>+(/- fia) ^ia]) />ia(r), (2) 

where Pio-(i*) stands for the orbital density, nicr(r) = 
pi(j{v)/ J dr^ Pier (1*0 denotes the normalized density, 
and En represent the exchange-correlation and Hartree 
energy functional, and v^c,a is the exchange-correlation 
potential. 

The pictorial interpretation of the non-Koopmans 
energy is simple. Graphically, the unrelaxed non- 
Koopmans energy Il^^{f) corresponds to the error made 
in approximating the ionization curve connecting the ini- 
tial state (/) to the final state (F) by a straight line 
whose slope is the derivative of the total energy at the 
transition point (Tf) (Fig. [3]). In fact, if the adiabatic 
Koopmans' theorem were satisfied, that is, if the energy 
€f^{f) were constant and equal to the ASCF energy dif- 
ference Ef^{fi(j) — £^^^(0), the ionization curve would be 
exactly linear and the non-Koopmans deviation Il^^{f) 
would vanish for any occupation / between and fi^. 

With the non-Koopmans measure in hand, it becomes 
possible to impose Koopmans' theorem, as shown in 
Ref. 0. Here, for the purpose of outlining the method, 
the non-Koopmans-corrected (NKC) functional can be 
described simply as a sum of weighted orbital constraints 
that are added to the Kohn-Sham total energy functional 
to penalize deviations from Koopmans' theorem: 

^NKC^^KS + ^^NKnu.KS(^NK)^ (3) 

ia 

In Eq. (|3]), the terms H^^^^ are the adiabatic non- 
Koopmans deviations that correspond to the chosen 
Kohn-Sham approximation, the coefficients afj^ repre- 
sent the weights of the non-Koopmans penalties, which 
must be calculated self-consistently to equate orbital en- 
ergies with ASCF energy differences [condition (a)] and 
the coefficients stand for the reference occupations, 
which must be fixed to preserve the accuracy of ASCF 
energies [condition (b)], as explained further below. 

At this point, it is important to note that Kohn-Sham 
density-functional theory is a ground-state theory, which 
is originally not intended to predict excited-state ener- 
gies. As a consequence, only the non-Koopmans energies 
of the highest occupied and lowest unoccupied molecu- 
lar orbitals can be rigorously defined within DFT since 
they correspond to ionization processes that take the N- 
electron ground state into the ground states with N — 1 
and A/" -hi electrons, respectively. In contrast, the ioniza- 
tion of other electron orbitals ends up into non-Aufbau 
states whose non-Koopmans energies are in principle not 
defined within DFT. 

As a consequence, the OD-DFT penalty sum that 
appears in Eq. ([3|) should be rigorously restricted to 
the frontier highest occupied and lowest unoccupied 
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states. Despite this recognized limitation, it is real- 
ized in practical computations that ASCF based upon 
local and semilocal DFT functionals predict accurate ex- 
citation energies provided that orbitals do not delocalize 
nonphysically.^^'^^ These computational observations in- 
dicate that generalizing non-Koopmans penalties to the 
full orbital spectrum in Eq. (jSj should yield a Koopmans- 
compliant OD-DFT method able to predict electronic 
structures with good accuracy. In practice, it is found 
that this generalization provides spectral predictions in 
very good agreement with experimental dataJ^^^^ 
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FIG. 4: Occupation dependence of the non-Koopmans- 
corrected total energy upon ionization of the 2p state of car- 
bon as a function of the penalty coefficient in the (a) un- 
relaxed and (b) relaxed cases. The underlying Kohn-Sham ap- 
proximation is LSDA. The reference occupation is of /IJf^ = |. 



Af/(an)-ANK(i)' 



(4) 



where A^^^(a) denotes the difference of orbital energies 
from (/) to (F) for the penalty coefficient afj^ = a. At 
that precise value, the energy curve is closely linear and 
Koopmans' theorem is satisfied. 

To illustrate this fact, we depict the influence of afj^ 
on the linearity of the energy curve upon ionization of the 
2p state of carbon in Fig. HJ Considering first the adia- 
batic case [Fig. ID^a)], it is seen that the linearity of the 
ionization energy, that is, the fulfilment of the adiabatic 
Koopmans' theorem, is achieved for 0^2^^^^ = 1. This result 
remains valid for any adiabatic process, as demonstrated 
by analyzing non-Koopmans residual errors. 

Instead, in the diabatic case, the value of afj^ that 
satisfies the generalized Koopmans' theorem is found to 
be lower than 1. As a matter of fact, for relaxed ioniza- 
tion of the 2p state of carbon, it is seen in Fig. HJ^b) that 
linearity is achieved for 0^2^^^^ close to |. Precisely, the 
converged value of the penalty coefficient that is deter- 
mined after 2 recursion steps [Eq. dU] is of 0.85. More 
generally, it has been shown that afj^ can be regarded 
as a relaxation coefficient that measures the stability of 
the final ionized state. Thus, for systems that consist of 
a lone electron around a closed shell (e.g., alkali- metal 
atoms), the outer-shell penalty coefficient equals 1. In 
contrast, afj^ becomes lower than 1 upon withdrawing 
an electron from a filled or partly filled shell (cf. Fig. 9 
in Ref. 28). 




To complete the presentation of the Koopmans- 
compliant functional, it remains to explain the deter- 
mination of the penalty coefficients afj^ and reference 
occupations f^^. The analytical determination of afj^ 
has already been presented. In this analysis, it has 
been shown that, starting from convex-energy local and 
semilocal DFT functionals, there always exist a value of 
afj^ lying between and 1 that equalizes the slope of 
the ionization curve (i.e., the orbital energies) at the ini- 
tial point (/) and that at the final point {F). Further- 
more, the value of afj^ that fulfills this condition can be 



FIG. 5: Occupation dependence of the non-Koopmans- 
corrected total energy upon relaxed ionization of the 2p state 
of carbon as a function of the reference occupation . The 
underlying Kohn-Sham approximation is LSDA. The penalty 
coefficient is set to be = 0.85, as calculated recursively 
from Eq. (g)). 

After explaining the calculation of the penalty coeffi- 
cients, we now turn our attention to the determination of 
the reference occupation. Based upon Slater's theorem 
and the sum rule satisfied by the exchange-correlation 
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hole, it can be shown that ff^^ should be set to ^ in 
order to preserve the accuracy of the underlying Kohn- 
Sham functional in predicting total-energy differences 
This fact is illustrated in Fig. [5] where it can be seen 
that the Koopmans-compliant energy difference matches 
the Kohn-Sham ASCF energy difference at half reference 
occupation for diabatic ionization of the 2p state of car- 
bon. Therefore, at the reference occupation /^^^ = 0.5 
and with the relaxation parameter 0^2^^^^ =0.85 that has 
been determined from Eq. (|4]), the Koopmans-compliant 
functional yields an orbital energy = 11.83 eV on a 
par with the accuracy of the ASCF prediction (11.72 eV) 
and experimental IP (11.26 eV). 

These results demonstrate the possibility of rectify- 
ing the nonlinear behavior of the ionization curve and 
restoring Koopmans' theorem [condition (a)] while pre- 
serving the accuracy of total-energy predictions [condi- 
tion (b)] within Koopmans-compliant OD-DFT. In ad- 
dition, Koopmans compliance presents the central ad- 
vantage of avoiding to treat negatively charge molecu- 
lar states [condition (c)] for both the evaluation of non- 
Koopmans penalties and the calculation of orbital levels, 
thereby allowing us to determine accurate acceptor levels 
straightforwardly, as shown in Sec. IIIIBi 

Further details on the implementation of one 
Koopmans-compliant method (the aNKCo computa- 
tional approach) are presented in the Appendix. 



III. RESULTS 

We now assess the predictive ability of Koopmans- 
compliant functionals in describing the donor-acceptor 
properties of a range of organic molecules relevant to 
OPV molecular junctions, first focusing on donor levels 
in Sec, nil Al and then examining the problem of acceptor 
levels in Sec. IIII B[ 



Donor levels 
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fluorene 





PTCDA phthalocyanine porphine 




tetraphenylporphyrin thiadiazole thiophene 



FIG. 6: Representative sample of organic molecules fRef. ITtI) 
for benchmarking the performance of electronic-structure 
methods. 



LSDA 
aNKCo 
Expt. 




FIG. 7: LSDA and aNKCo donor levels compared with ex- 
periment (cf. references in Table [l]) for the organic molecules 
listed in Fig. [6l 



The initial benchmark that is depicted in Fig. [6] con- 
sists of a representative set of molecules, which has been 
studied in Ref. 117|. These molecules constitute elemental 
components for relevant photovoltaic organics and allow 
us to assess the accuracy of the self-interaction method 
in capturing the electronic structure of organic materials 
across a representative sample of chemical compositions 
and molecular sizes. 

Koopmans-compliant calculations are carried out at 
the aNKCo level. We employ conventional norm- 
conserving pseudopotentials to represent atomic cores 
with a cutoff energy of 60 Ry in expanding wave func- 
tions. All of our calculations take into account spin po- 
larization. Employing auxiliary-function correction^ a 
vacuum separation of 9 A is sufficient to achieve conver- 
gence of the electronic structure within 0.15 eV, even for 
the large PTCDA and phthalocyanine molecules. In all 



cases, less than two recursions of Eq. (j4]) are needed in 
converging the penalty coefficient a of the highest occu- 
pied state of the molecule, and the same penalty weight 
is imposed on the full spectrum (as justified by the sen- 
sitivity analysis presented in the Appendix). Damped 
electronic dynamics is used in optimizing electronic de- 
grees of freedom with the inner-loop procedure described 
in Ref. 60 and with a convergence threshold of 10 ~^ Ha 
on the electron-dynamics kinetic energy. Regarding the 
experimental literature, it is important to mention that 
all of the measurements that are reported here corre- 
spond to molecules isolated in the gas phase. Therefore, 
environmental effects arising from a surrounding solvent 
or from an embedding bulk need not be included in our 
calculations. 

LSDA and aNKCo predictions for the frontier donor 
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TABLE I: Absolute donor levels of benchmark organic 
molecules within LSDA, PZ, aNKCo, and GW compared with 
experimental vertical ionization energies (cf. Ref. 45i and ref- 
erences below). 

Donor levels (eV) 





LSDA 


PZ 


aNKCo 


GW" 


Expt. 


benzothiadiazole 


6.40 


10.29 


9.03 


8.56 


9.15^ 


benzothiazole 


6.20 


9.66 


8.64 


8.48 


8.85" 


fluorene 


5.64 


8.69 


7.51 


7.64 


7.91^ 


PTCDA 


6.32 


9.66 


8.12 


7.68 


8.2" 


phthalocyanine 


5.21 


8.31 


6.54 


6.10 


6.41^ 


porphine 


5.24 


8.34 


6.66 


6.70 


6.9^ 


tetraphenylporphyrin 


4.89 


8.10 


6.21 


6.20 


6.39^ 


thiadiazole 


7.13 


11.73 


10.10 


9.89 


lo.ir 


thiophene 


6.02 


10.97 


8.76 


8.63 


8.85^' 


MAD 


2.20 


1.43 


0.16 


0.32 




RMS 


0.64 


0.45 


0.11 


0.15 





"Reference 
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Reference 46. 
"Reference 47 
"^Reference 48 
"Reference 49. 
•^Reference 5Q. 
^Reference [51 
^Reference [5^ 
^Reference 53. 
^Reference 54. 
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FIG. 8: LSDA and aNKCo donor levels compared with ex- 
periment for benzene (Ref. ..55., ) and oligoacenes (Ref. [56l ). 



levels (i.e., the absolute energy of the highest occupied 
orbital) of benchmark molecules are reported and com- 
pared with experiment in Fig. [71 (Here and throughout 
the comparative assessment, we focus on LSDA results 
rather than on GGA predictions due to the fact that 
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FIG. 9: LSDA and aNKCo donor levels compared with ex- 
periment for fullerenes. Experimental vertical ionization po- 
tential are from Ref. IstI for C50, Ceo, and Cso (Dsd) and from 
Ref. m for C70. 



aNKCo is applied to LSDA; this permits to directly re- 
solve the effective influence of aNKCo on the underlying 
functional. Note that GGA orbital energies would be 
very close to the LSDA results on the scale of the self- 
interaction error.) 

First, we observe a marked underestimation of absolute 
donor levels within LSDA. The underestimation is above 
1 eV for porphyrin complexes and reaches 2.5-3 eV for 
strong aromatic donors. This incorrect trend is rectified 
by aNKCo that reduces the error down to a few tenths 
of an eV. Furthermore, the accuracy of aNKCo compares 
favorably with that of the original self-interaction cor- 
rection (PZ)^^ and GW approximation. Indeed, data 
reported in Table |I] reveal that PZ provides some im- 
provement over LSDA, reducing the mean absolute error 
from 2.20 eV to 1.43 eV. Much better agreement with 
experiment is achieved within self-consistent GW^^ with 
a mean absolute error as low as 0.32 eV and a stan- 
dard deviation of the error of 0.15 eV. Despite the very 
good performance of GW calculations,^^ the Koopmans- 
compliant aNKCo method is found here to be more pre- 
cise in predicting the highest donor levels with an error 
of 0.16 eV and a low standard deviation of 0.11 eV. 

After highlighting the remarkable precision of aNKCo 
across a variety of chemical compositions, we now assess 
its ability to capture the influence of molecular size on the 
electronic structure. For this complementary benchmark, 
we consider two important families of organic molecules, 
namely, acene oligomers and fullerene clusters. We em- 
ploy the same calculation parameters as those above — 
with the exception of the plane-wave cutoff that can be 
reduced to 40 Ry without altering numerical convergence. 

For oligoacenes (Fig. [8]), we observe that LSDA is in er- 
ror of 2.8 eV for the donor level of benzene. The error de- 
creases gradually with the length of the chain but is still 
as large as 2 eV for hexacene. In terms of relative errors. 
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TABLE II: Absolute acceptor levels of benchmark organic 
molecules within LSDA, PZ, aNKCo, and GW compared with 
available experimental data. 



TABLE III: Absolute acceptor levels of acene chains of in- 
creasing length within LSDA, aNKCo,and GW compared 
with experimental data. 



Acceptor levels (eV) 



Acceptor levels (eV) 





LSDA 


PZ 


aNKCo 


GW" 


Expt. 






LSDA 


aNKCo 


GW" 


Expt. 


benzothiadiazole 


3.52 


3.42 


1.07 


0.42 




benzene 




1.38 


0.03 






benzothiazole 


2.35 


2.20 


<0 


<0 




naphthalene 


2.27 


0.49 






fluorene 


2.05 


2.00 


<0 


<0 




anthracene 




2.85 


1.16 


0.29 


0.53^ 


PTCDA 


4.80 


4.81 


3.19 


2.68 




tetracene 




3.22 


1.60 


0.93 


1.07" 


phthalocyanine 


3.79 


3.82 


2.34 


2.07 




pentacene 




3.48 


1.98 


1.36 


1.39" 


porphine 


3.28 


3.19 


1.53 


1.39 




hexacene 




3.68 


2.18 






tetraphenylporphyrin 


3.07 


3.04 


1.68 


1.49 


1.69(10)^ 


octacene 




4.27 


2.48 






thiadiazole 


2.95 


2.78 


<0 


<0 
















thiophene 


L59 


1.38 


0.04 


<0 




MAD 




2.19 


0.59 


0.13 
















RMS 




0.09 


0.04 


0.09 




^Reference fl^. 
^Reference \6u. 












"^Reference 
^Reference 
^Reference 


17 
62 
63 





the LSDA underestimation fluctuates in the narrow range 
29-31%. Applying aNKCo restores the agreement be- 
tween electronic- structure calculations and experiment, 
reducing the error to less than 0.2 eV regardless of the 
length of the chain. Similarly, for fullerenes (Fig. [9]), 
aNKCo brings about some considerable improvement in 
capturing the delicate influence of cluster size. Indeed, 
despite the fact that the LSDA error varies in a much 
wider range than for acene oligomers (22-28%) and that 
the dependence of the highest occupied level as a func- 
tion of the molecular size is clearly nonmonotonic, the 
QfNKo error never exceeds 0.2 eV compared to available 
experimental data, providing further conflrmation of the 
predictive precision of the Koopmans-compliant method 
in describing subtle electronic-structure trends. 

In addition, an important advantage of the aNKCo 
method lies in its moderate computational cost. In fact, 
converging the full electronic spectrum of an extended 
PTCDA molecule with a large plane-wave cutoff of 60 
Ry on a conventional 16-processor machine requires half 
of a CPU day, whereas the same calculation within HF 
and hybrid-DFT is 3 times more costly using compa- 
rable CP implementations. Indeed, as explained previ- 
ously, although conventional self-interaction corrections 
and hybrid-DFT both scale quartically [0{N^)] with the 
size of the electronic system, the algorithmic prefactor 
of self-interaction methods is lower than that of hybrid- 
DFT since the former do not require to evaluate pairwise 
exchange contributions. 



B. Acceptor levels 

After validating the aNKCo method for occupied 
donor levels, we now concentrate on the prediction of un- 
occupied acceptor levels. As already mentioned, accep- 



tor levels represent a central difficulty for ASCF predic- 
tions due to the inability of conventional DFT methods 
to properly stabilize negatively charged states. In what 
follows, we demonstrate that aNKCo provides a solution 
to this important methodological limitation. 

In calculating acceptor levels, we use the same com- 
putational parameters as those of donor-level calcula- 
tions. In particular, the aNKo penalty weights remain 
unchanged. First-principles predictions for benchmark 
organic compounds are reported in Table |TT1 Our results 
clearly confirm the poor performance of LSDA and PZ 
in predicting acceptor levels; both methods overstabilize 
acceptor states by several eVs relative to GW due to self- 
attraction. In contrast, aNKCo and GW predictions are 
found to be in close agreement. In quantitative terms, the 
discrepancy between aNKCo and GW is as low as 0.15 
eV for porphine and is at most of 0.6-0.7 eV for benzoth- 
iadiazole. It should also be noted that our comparison 
reveals that aNKCo acceptor levels are always more sta- 
ble than their GW counterparts. Additionally, although 
experimental data for benchmark molecules are scarce, 
comparison with the measured vertical electron affinity 
of tetraphenylporphyrin suggests that the predictive ac- 
curacy of aNKCo is comparable with the accuracy of GW 
predictions. 

To gain further insight into the performance of first- 
principles methods in predicting acceptor levels, we con- 
centrate on the lowest unoccupied state for acenes (Ta- 
ble llllj) and for fullerenes (Table ITV|) . For weak acene 
acceptors, the first notable feature is the propensity of 
LSDA to overestimate absolute acceptor levels. This ex- 
pected nonphysical tendency is rectified to a large ex- 
tent by aNKCo, which reduces the mean absolute error 
from 2.19 to 0.59 eV relative to vertical ionization exper- 
iments. Nevertheless, it is seen that GW performs sig- 
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TABLE IV: Absolute acceptor levels of fullerene clusters of in- 
creasing size within LSDA, aNKCo, QMC (quantum Monte- 
Carlo), GWo (nonself-consistent), scGWf (self-consistent 
with vertex corrections) compared with experimental data. 



Acceptor levels (eV) 





LSDA 


aNKCo 


QMC" 


GWo" 


scGWf" 


Expt. 


C20 


4.24 


2.19 


1.76(11) 


3.55 


2.36 


2.25^ 


C24 


5.06 


2.88 


2.57(11) 


4.19 


2.88 




C50 


5.20 


3.49 


3.52(14) 


4.75 


3.73 


3.10' 


Ceo 


4.27 


2.64 


2.23(19) 


3.87 


2.98 


2.69^ 


C70 


4.10 


2.61 


2.46(11) 


3.98 


2.83 


2.76" 


Cso (Dsd) 


5.00 


3.91 


3.25(10) 


4.62 


3.88 


3.70' 


Cso (Ih) 


4.56 


2.99 


3.90(11) 


5.17 


4.38 




MAD 


L66 


0.18 


0.42 


1.25 


0.26 




RMS 


0.32 


0.12 


0.06 


0.26 


0.21 





"Reference 
^Reference 65. 
'Reference 66.. 
^Reference [63. 
'Reference ^6§. 



nificantly better than aNKCo, bringing the error down 
to 0.13 eV. In the same vein, we note that the naph- 
thalene acceptor level is predicted to be stable within 
QfNKCo in contradiction to experiment. In fact, among 
all of our benchmark calculations, acenes represent the 
only instance where GW performs better than aNKCo, 
suggesting that GW should be preferred over aNKGo in 
capturing shallow acceptor levels. 



IV. CONCLUSION 

In order to overcome the predictive deficiency of con- 
ventional DFT methods in capturing the donor and ac- 
ceptor levels of molecular complexes, we have presented 
a correction procedure that enables one to eliminate the 
nonphysical curvature of the total energy upon with- 
drawing (injecting) electrons from (into) the molecular 
system. When applied to LSDA, the procedure yields 
an orbital-dependent functional that fulfills the gener- 
alized Koopmans' theorem, thereby restoring or impos- 
ing the agreement between frontier orbital levels and 
ASCF energy differences, and overcoming central limi- 
tations in applying ASCF techniques to predict accep- 
tor levels. Specifically, the aNKGo formulation of the 
Koopmans-compliant method has been applied to OPV 
compounds spanning a wide range of chemical composi- 
tions and molecular sizes, thereby highlighting the con- 
ceptual simplicity and computational performance of the 
method with an accuracy comparable to that of many- 
body perturbation theory. 

Although the study has essentially focused on de- 
termining frontier donor and acceptor levels, the non- 
Koopmans corrective procedure can be straightfor- 
wardly extended to compute the full electronic spec- 
trum of molecular systems in a single first-principles 
calculation.^^ For future studies, we plan to apply the 
Koopmans-compliant method to other relevant fami- 
lies of photovoltaic compounds, including metallic por- 
phyrins and phthalocyanines, and to generalize the ap- 
proach to describe strongly interacting dyads and bulk 
organic semiconducting materials at both the electronic 
and excitonic levels. 
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Appendix A: Computational details 

In this appendix, we present the computational im- 
plementation of the Koopmans-compliant method in the 
CP (Car-Parrinello) code^^ of the quantum-espresso 
open-source project. Specifically, we describe the cal- 
culation of contributions to the total energy and effective 
Hamiltonian that arise from the non-Koopmans correc- 
tion, we explain algorithmic choices in optimizing the 
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electronic structure, and we detail the evaluation of the 
penalty coefficients. 

Before proceeding to the calculation of non-Koopmans 
contributions, it should be noted that CP is a plane- wave 
code. Adopting a plane- wave representation here allows 
us to envision applications to relevant periodic molec- 
ular structures, such as semiconducting polymers and 
bulk materials. Despite the benefit of the plane- wave 
approach, the evaluation of electrostatic contributions 
to the non-Koopmans energies [Eq. (|2])] in the context 
of plane waves requires to correct periodic-image errors 
that arise from the use of the supercell approximation 
in describing nonperiodic systems. To eliminate arti- 
ficial periodic-image interactions, we employ auxiliary- 
function corrections that consist in removing singular- 
ities in the reciprocal-space summation of electrostatic 
terms through the addition of a regularizing auxiliary 
functions to the point-charge electrostatic kernel. The 
choice of this reciprocal-space scheme is motivated by 
the fact that it allows us to use fast-Fourier-transform 
(FFT) techniques, thereby minimizing the computational 
burden associated with the repeated calculations of the 
electrostatic energy and potential for each orbital of the 
system; without this correction the calculation would be 
prohibitively expensive. 

The second difficulty in calculating non-Koopmans 
terms lies in the assessment of the second-order func- 
tional derivatives that appear in the expression of 
the variational contributions to the effective potential, 
namely, the and w^^^^ terms defined in Ref. i28i. 

To address the calculation of these terms, we employ the 
computational subroutines written by Dal Corso and de 
Gironcoli in the context of the DFPT (density-functional 
perturbation theory) calculation of phonon dispersions.^^ 
These subroutines are based upon explicit expressions 
of second-order derivative contributions and also em- 
ploy finite-difference techniques. Similar subroutines dis- 
tributed as independent libraries (e.g., the LIBXC module 
of the OCTOPUS code) ''^^''^ can be conveniently ported to 
compute second-order terms. 

Overall, the calculation of non-Koopmans contribu- 
tions to the total energy and effective Hamiltonian is 
already less expensive than that of exchange contribu- 
tions in hybrid-DFT functionals.^^ Nevertheless, in order 
to gain further computational performance, it is neces- 
sary to reduce the computational burden associated to 
the calculation of and w^^-^. To this end, a suc- 

cessful strategy consists in not updating the reference 
density pfj^{r) = p{r) + (/j^^ - fia)nia{r) that appears 
in the expression of non-Koopmans energies at every iter- 
ation of the self-consistent electronic-structure optimiza- 
tion but only periodically for each given number of itera- 
tions, thereby avoiding to calculate and '^^^^^r 
each self-consistent step. This computational approach 
that consists in neglecting the and '^^d contribu- 
tions to the self-consistent OD-DFT potential is referred 
to as the NKCq method. 

Remarkably, it is found that relaxing self-consistency 



via the NKCq method preserves and actually improves 
the accuracy of predicted electronic properties relative 
to NKC for difficult molecular complexes due to its in- 
creased localization strength, with the caveat that only 
orbital properties and not total energies can be evalu- 
ated using NKCq.^^ Variational alternatives to NKCq, 
endowed with a comparable localization strength can be 
derived at the cost of increased sophistication. A pre- 
liminary discussion of variational non-Koopmans self- 
interaction corrections is presented in Ref. i7^. In the 
present study, we focus on NKCq due to its simplicity 
and accuracy in describing individual organic donors and 
acceptors. Admittedly, more work would be necessary for 
the variational description of systems where strong elec- 
tronic delocalization prevails. 

Regarding computational implementation, NKCq is 
straightforward to program since it simply relies on sub- 
stituting the DFT Kohn-Sham potential with the OD- 
DFT potential written below: 

i (Al) 

(l-aNt<KSj^]+«NK,KSj^NK] 

where v^^^ denotes the sum of the Hartree and exchange- 
correlation potentials and pfj^ has been defined above. 
Once the subroutines that calculate v^^^ are identified, 
this modification can be performed with some minimal 
coding experience. 

We now turn our attention to the optimization of the 
electronic structure. In this regard, we first note that 
only few of the conventional minimization algorithms ap- 
plicable to DFT can be employed in the context of the 
OD-DFT corrections. In particular, conventional itera- 
tive diagonalization would represent a very poor choice 
since eigenfunctions of different Hamiltonians would have 
to be calculated for each of the orbitals in the system, 
thereby considerably increasing computational complex- 
ity. Instead, direct-minimization approaches, such as 
conjugate-gradient or damped-dynamics techniques, are 
more suited to OD-DFT calculations. Additionally, since 
OD-DFT self-interaction corrections lead to energy func- 
tional that are variant with respect to unitary rotation 
of the orbital manifold contrary to DFT, it is also nec- 
essary to consider unitary-rotation degrees of freedom in 
minimizing the energy. Several approaches have been 
proposed and can be used to this end.^^^^^^^ 

The final hurdle in implementing the Koopmans- 
compliant method pertains to evaluating the penalty co- 
efficients afj^. 

Here, one first source of computational savings springs 
from the fact that afj^ is much less sensitive to the 
resolution of the plane-wave grid than orbital-level pre- 
dictions. In quantitative terms, in the case of fur an 
(Fig. [To]), reducing the cutoff from 60 Ry to 15 Ry causes 
the relaxation coefficient to vary by 3%, whereas such a 
low plane- wave cutoff in expanding wave functions would 
be intolerably low for the description of orbital levels. 
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Therefore, since afj^ can be evaluated on coarse plane- 
wave grids, the preliminary determination of the weights 
of the non-Koopmans penalty afj^ represents a marginal 
fraction of the computational cost associated with the 
calculation of orbital levels. 




10 20 30 40 50 60 

wave-function cutoff energy Ecut (Ry) 

FIG. 10: Non-Koopmans penalty coefficient for the highest 
occupied molecular orbital [Ia2(7r3)] of furan as a function of 
wave-function cutoff energy. 

The second important source of cost reduction in prac- 
tical calculations is the fact that the relaxation coeffi- 
cients afj^ varies in a limited range, typically, between 
0.5 and 1. To illustrate this fact, we report the value of 
the coefficients afj^ for the electronic states of carbon 
in Fig. [TTJ On the basis of these calculations, we can 
assess the sensitivity of electronic-level predictions as a 
function of the afj^^s. Indeed, LSDA is typically in error 
of A = 40% in predicting the electronic levels of atoms 
and molecules. As a consequence, a simple sensitivity 
estimate reveals that an error of a — afj^ in the determi- 
nation of the self-consistent penalty coefficient translates 

into an error of — in the electronic level relative 

to the exact Koopmans-compliant estimate (that is, rel- 
ative to the ASCF energy difference). Thus, in the case 
of carbon, making the approximation of equating all of 
the penalty coefficients to that of the highest occupied 2p 
state (i.e., a = 0.85), we obtain errors as low as 2% for 
the spin-up 2s state and lower than 14% for the spin-down 
Is state, which in both cases represents a significant re- 
duction of the LSDA error relative to ASCF predictions. 
This sensitivity analysis provides a clear justification 



for the QfNKCo method that consists in setting the 
penalty coefficients afj^ to be all equal to the same value, 
thereby avoiding in particular to treat negatively charged 
states. In explicit terms, aNKCo consists in performing 
the following substitution in lieu of Eq. ()Aip : 



(1 



«)«Hxc 



i 
[P] 



(A2) 



To complete the presentation of aNKCo, it is relevant 
to discuss the size consistency of the method. Within 
QfNKCo, size consistency may be lost when considering 
a molecular system composed of separate fragments; in 
this situation, the a coefficient calculated for the global 
system may differ from the a coefficients that correspond 
to each of its subparts. (The same problem is known to 
arise in range-separated hybrid-DFT approaches where 
a single optimally tuned parameter 7 is employed to im- 
pose Koopmans' theorem.) An immediate solution to ad- 
dress this problem is to resort to NKCo, which preserves 
the size consistency of the underlying local or semilocal 
functional owing to the orbital-specific coefficients dfj^- 
Note that other simple solutions may be employed, such 
as the approach that consists in using a different a for 
each separate molecular fragment. Following this proce- 
dure, orbit als localized on a given molecular fragment feel 
the same a while orbitals localized on different fragments 
have different a's, thereby preserving size consistency. 
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FIG. 11: Penalty coefficients afj^ that quantify orbital relax- 
ation upon ionization for the orbitals of carbon in the spin-up 
and spin-down channels. 
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